/*!
 * \file laplace1D.cpp
 * Solving the Laplace equation in 1D.
 ***************************************************************************
 *  \author    Luis M. de la Cruz Salas [ Fri Mar 21 12:26:33 GMT 2010 ]
 ***************************************************************************
 * This code solves the Laplace equation in one dimension:
 * \f$ \nabla^2 T = 0 \f$.
 * This equation is solved in the interval \f$ x \in [0,1] \f$.
 * The initial condition for the temperature is \f$ T = 10 \cos(3 \pi x) \f$
 * and the boundary conditions: \f$ T=10 \f$ for \f$ x=0 \f$, and \f$ T=-10 \f$ 
 * for \f$ x = 1 \f$.
 * You can choose between Gauss-Seidel or Thomas algorithms to solve this 
 * problem. In the case of Gauss-Seidel you need to specify the maximum 
 * number of iterations.
 *  \par Compiling and running
 *  Modify the variables BASE and BLITZ in the file \c tuna-cfd-rules.in 
 *  according to your installation and then type the next commands:
 *  \verbatim 
   % make
   % ./laplace1D  \endverbatim
 * 
 * The results are stored in \c Data1D.
 ***************************************************************************
 */
#include <iostream>
#include <cmath>
#include "Meshes/Uniform.hpp"
#include "Equations/ScalarEquation.hpp"
#include "Schemes/LaplaceScheme.hpp"
#include "Solvers/TDMA.hpp"
#include "Utils/num_utils.hpp"
#include "Utils/inout.hpp"
#include "Utils/GNUplot.hpp"
using namespace Tuna;

// Scalar field in one dimension
typedef TunaArray<float,1>::huge ScalarField1D;

timer crono;

int main()
{
    float length_x = 1.0; 
    int num_nodes_x, max_iter = 0;
    bool solver;
    cout << "\n No. nodes = "; cin >> num_nodes_x;
    cout << "\n Solver : \n   (0) Gauss-Seidel \n   (1) Thomas \n option : ";
    cin >> solver;
    if (!solver) {
      cout << "\n Max_iter Gauss-Seidel = "; cin >> max_iter;
    }

    StructuredMesh<Uniform<float, 1> > mesh(length_x, num_nodes_x);
    mesh.print();
    float dx = mesh.getDelta(X);

    ScalarField1D T (mesh.getExtentVolumes());
    ScalarField1D Tn(mesh.getExtentNodes());

    firstIndex i;
    Tn = 10 * cos (3 * PI * dx * i);
    NumUtils::interpolateToVolumes(T, Tn);
    InOut::writeToFile(Tn, 0, "./Data1D/temp.", dx);

#ifdef WITH_GNUPLOT
    GNUplot plotter;
    plotter("set grid");
    plotter("p \"./Data1D/temp.00000\" title \"Initial T\" w lp");
    std::cout << "\n >---> Press <1> and then <enter> to continue \n";
    int pause; std::cin >> pause;
#endif

    SparseMatrix< Diagonal<float, 1> > A(num_nodes_x); 

    ScalarField1D                b(num_nodes_x);

    ScalarEquation < LaplaceScheme<float, 1> > laplace(T, A, b, 
						       mesh.getDeltas());
    laplace.setGamma(1.0);
    laplace.setDirichlet(LEFT_WALL);
    laplace.setDirichlet(RIGHT_WALL);
    laplace.print();

    laplace.calcCoefficients();

    if (solver) {
      Solver::TDMA1D(laplace);
    } else {
      cout << "\n | N. iter = " << Solver::GaussSeidel1D(laplace, 0.001, max_iter);
    }
    laplace.update();

    NumUtils::interpolateToNodes(Tn, T);
    InOut::writeToFile(Tn, 1, "./Data1D/temp.", dx);

#ifdef WITH_GNUPLOT
    plotter("set grid");
    plotter("p \"./Data1D/temp.00000\" title \"Initial T\" w lp, \"./Data1D/temp.00001\" title \"Final Sol\" w lp");
    std::cout << "\n >---> Press <1> and then <enter> to continue \n";
    std::cin >> pause;
#endif

    return 0;
}

